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Abstract 

A method of deriving quadrature rules has been developed which gives 
nodes and weights for a Gaussian-type rule which integrates functions of 
the form: 

ff .\ aix,y,t) b{x,y,t) ,\i w +^2 , 21I/2 , 

" {x ~ ty + y2 + [(^ _ tf + y2] V2 +^(^' ^) l»g[(^-^) +2^ 1 -^d{x,y,t), 

without having to explicitly analyze the singularities of /(a;, j/, t) or sepa- 
rate it into its components. The method extends previous work on a sim- 
ilar technique for the evaluation of Cauchy principal value or Hadamard 
finite part integrals, in the case when y = 0. The method is tested by eval- 
uating standard reference integrals and its error is found to be comparable 
to machine precision in the best case. 

1 Introduction 

An important part of the application of the boundary element method (BEM) to 
physical problems is the calculation of the field, from the solution on the bound- 
ary. This requires the evaluation of integrals which contain 'almost-singular' 
integrands which are not properly handled by standard Gaussian quadratures. 
For example, if we consider a two-dimensional potential problem, such as the 
Laplace or Helmholtz equations: 

V^^ = 0, (la) 
V20-Hfc^0 =0, (lb) 

where fc is the wavenumber, the potential at some point in the field will be: 
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Figure 1: Geometry of boundary element: collocation points indicated by filled 
circles, field point position x = Rcos9, y = Rsm9. 

where G(x, xi) is the Green's function for the problem, F is the boundary of the 
domain with normal n and variables of integration are indicated by subscript 
1. In the case of the Helmholtz equation, the singular behaviour of the Green's 
function will be related to the Green's function of the corresponding Laplace 
equation [8, for example]. The Green's functions for the Laplace equation will 
have a logarithmic singularity in the planar and axisymmetric case (where the 
Green's function is proportional to an elliptic integral) and also in the case of 
an asymmetric problem in an axisymmetric domain [3,5]. 

The Green's function for the two-dimensional Laplace equation is G = 
log |x — xi]. If the boundary integral problem has been solved using the BEM, 
the evaluation of the potential in the field gives rise, on each element, to an 
integral of the form: 

/W(x,y)= j'^{\og\{x-tf+y^\+f{x,y,t))L{t)dt, (3) 

where t is the local coordinate on the element, L{t) is a shape fmiction, typically 
a polynomial, and x and y are the field point coordinates in the local coordinate 
system, figure 1. This integral is non-singular but suffers from an 'offstage 
singularity' [1] if the field point is near the element, so that the argument of the 
logarithm becomes small. 

If it is required to determine the first or second derivatives of the potential, 
for example, to determine a velocity or velocity gradient in a fluid-dynamical 
problem, integrals of the form: 

I^'Hx,y) =/^^p-^m_^+i(f)g(^,y,t)df, (4a) 

I^^\x,y) =J^^^J^+L{t)h{x,y,t)dt, (4b) 

arise. By analogy with the case when y = 0,we refer to these integrals as 'near- 
singular' and 'near-hypersingular' respectively. When y = 0, the field point 
lies on the boundary and the integrals must be treated as Cauchy principal 
values [12, page 37] or as a Hadamard finite-part [12, page 31]. In this C&SC^ cl 
number of procedures exist for the accurate evaluation of the integrals [4, 13, 14, 
for example] 

A recently-developed method [6, 11] develops quadrature rules which can be 
used as 'plug-in' replacements for standard Gaussian quadratures in cases where 
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the integrand contains singularities of the form of equations (3) and (4). An 
important feature of these rules is that they require no analysis of the integrand 
before use. When the potential is differentiated, yielding equation (4a), the 
function g{x,y,t) will usually contain logarithmic singularities which must be 
properly accounted for, as well as the leading singularity. In many problems, 
the Green's function will be complicated and the explicit identification of the 
singularities will be time-consuming. It is preferable to have a qiiadrature rule 
which correctly integrates all of the singularities present without requiring a 
detailed analysis. The previously published method [6, 11], generates such rules 
for singular and hypersingular integrals. This paper extends this method to 
'near-singular' and 'near-hypersingular' integrals, giving a simple, easily imple- 
mented technique producing quadrature rules which are direct replacements for 
Gaussian quadrature in BEM codes. 



2 Quadrature rules for near-singular integrals 

Before developing the quadrature rules needed for near-singular integrals, it is 
worth examining the reason for the breakdown of standard Gaussian quadrature. 
This will also give an indication of when specialized rules should not be used 
and standard quadratures are better. Consider the integral: 

I = /_! [(x-t)2+y2]i/2 di, (5) 
= i f 1 1 M 

RJ-1 ((t/fi)2-2tcos6l//<+l)i/2 ""-I 

where R'^ = + and 9 = ta,n~^y/x. This can be expanded using the 
generating function for Legendre polynomials [9, 8.921]: 

^ EZot'Pki^)^ |i| <nnn|z±(z2_l)V2|^ (6) 



(l-22i + t2)l/2 

EZot-^'^'^Pkiz), \t\>ma^\z±{z'-l)'/% (7) 

to yield: 

=EZo^PUcose), t<R (8a) 

= EZoWTTPk{cos0), t>R. (8b) 

An A^-point Gaussian quadrature rule integrates exactly polynomials up to order 
2N — 1. When R^ 1, the integrand is well approximated by equation (8a) and 
the estimate of / returned by Gauss-Legendre quadrature: 

^-2E^^^,^^^P..(cos^), (9) 

will be accurate. If, however, R < 1, part of the integrand will be given by 
the inverse power series of equation (8b) which cannot be correctly integrated 
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by the standard Gaussian quadrature. Likewise, even if i? > 1 but is not large 
enough to make the terms of equation (8a) decay fast enough, there will be a 
large error in the estimate of I. Looking ahead to the results presented in §3, it 
is expected that standard Gauss-Lcgcndre quadrature will give good rcsiilts for 
large R and/or in cases where P„(cos6') is small. Otherwise, a specialized rule 
will be necessary. 



2.1 Evaluation of quadrature rules 

The algorithm to be developed gives an iV-point rule which integrates a function 
of the form: 

/(-'^'*) = f ,2 + 1% t) log[(^-^)^ V] '/'+d{x, y, t), 

(10) 

where a, b, c and d are taken to be polynomials of order up to M and the 
integral 

/•I ^ 

where U are the quadrature points of an TV-point Gauss-Legendre quadrature 
and the weights Wi are to be determined. A previously-developed algorithm [6, 
11] for the computation of quadrature rules gives a method for the evaluation 
of Wi when y = 0. The approach is conceptually simple — the weights are found 
as the solution to the system of equations: 



.1 N 

/ f{x,y,t)dtfv'^Wif{x,y,ti), (11) 



Y,[jl)ij]wj = rrii, i=l,...,iM, (12) 



where tpij are the weighted Legendre polynomials at the quadrature points tj: 

P^-i{t,) l<i<M, 
F,_M-i(ij)log[(:r-t,)2+j/2]i/2 M+l<i<2M, 

P^-2M-l{t3)/[{x - tjf + y2]i/2 2M + I < I < 3M, 



Pi_^M-i{tj)[{x - tjf +2/2] 3M + 1 < i < 4M. 
and the moments mi are the integrals of V'i 

' /^^Pi_i(t)dt l<i<M, 

/^^ P^-M-l{t) l0g[(,T - t)2 + y2]l/2 M + 1<1< 2M, 

/^,P,_2M-i(i)/[(x-t)2+2/2]i/2dt 2M + l<t<3M, 
, J'^ Pi-3M-i{t) [{x -tf + y^]dt 3M+l<i< AM. 



mi= < 



(14) 



The system of equations is solved using the appropriate LAPACK solver [2], 
in the least squares sense when N > AM and in the minimum norm sense 
when A'' < 4M. The only issue which must be clarified is the evaluation of the 
moments, mj. In the case when y = 0, the integrals are true Cauchy principal 
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values or Hadamard finite parts and there exist formulae for their evaluation 
in terms of associated Legendre functions [10] or simple finite part integrals 
combined with exact Gaussian quadratures [6]. In this case, however, such 
simple formulae are not available and a different approach is required. 



2.2 Integration of weighted Legendre polynomials 

The method outlined in §2.1 for the evaluation of the quadrature weights Wi 
requires the evaluation of the moments rrii where 



mn = j u{t)Pn{t)dt, (15) 

where weighting function u{t) will be one of log[(x — t)"^ + y"^]^^"^ or [{x — t)^ + 
j^2j-ra/2^ n = 1 or 2. A general method can be applied to finding the integrals 
of weighted Legendre polynomials, using basic functional relations and simple 
formulae for the integrals of elementary functions, readily found in standard 
references [9]. 

Assuming that integrals of the form 



t"-u{t) dt, (16) 



can be evaluated (the formulae required for this paper are given in the appendix), 
the expansion of powers of t in terms of Legendre polynomials [9, 8.922.1]: 

= i^hPoit) + ELi(4fc + 1) ^^^r^^s;'::^':^^^, P2k{t) (m) 

can be used to show that: 

J2n+1 = 27^/^1 ^'l(OwWdi 

+ ELi(4fc + 3) (.iSpifeSl^ P2,MtMt) dt. (18b) 

Then, given the integrals J^, i = 0, 1, . . . , A^, the corresponding integrals of the 
weighted Legendre polynomials, m,;, can be evaluated via: 

Ci{2n+l)m2„ = J2„- Jo/(2n+l)-E"=i(4j + l)C,m2,., (19) 



mo -Jo, Cj - 2l+2j+3^j-'^' ^0-2^+3 



and 



A(2n + 3)m2„+i = J2„+i - 3Ji/(2n + 3) - E"=i(4j + 3)£»,m2,+i, (20) 
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The evaluation can be carried out in place, with values of nii overwriting J^. In 
double precision, the procedure gives accurate results for n ^ 32 before overflow 
errors cause problems. For BEM applications, where the shape functions are 
typically of order no greater than 3, this causes no special difficulties, but would 
limit use of the method in other areas. Obviously, for any application other than 
straight elements with G = log |x — xi |, a rule with M > 3 will be needed, but 
it is unlikely that M will need to be greater than 32. 

2.3 Summary of method 

To summarize, the procedure for computing a quadrature rule with N points 
which can integrate functions of the form of equation (10) where a, b, c and d 
are polynomials of order up to M is as follows: 

1. find the quadrature points ti for an A^-point Gaussian quadrature, using, 
for example, the method of Davis and Rabinowitz [7]; 

2. evaluate the weighted Legendre polynomials ipij at each of the quadrature 
points; 

3. compute the moments rrii using the method of §2.2; 

4. solve [ipij]wj = rrii using an appropriate solver. 
The integral of f{x,y,t) is then approximated by: 

.1 N 

/ f{x,y,t)dt^y2wa{x,y,ti). (21) 
J-^ i=0 

3 Numerical tests 

The quadrature method developed in §2 is tested by applying it to the evaluation 
of a reference integral. Before carrying out these tests, it is of some interest to 
examine the behaviour of the quadrature weights Wi with respect to the field 
point position. Figure 2 shows the root mean square difference 5 between the 
rule of this paper with A'' = 64, M = 16 and a standard 64-point Gaussian 
quadrature, with 

1/2 

(22) 

where Vi are the weights of the standard rule. 

The difference between the rules is shown as a function of R and 9. It is clear 
that close to the element {R and/or 9 small), the weights of the new quadrature 
are very large, as they have to cancel the large value of the integrand close 
to the near-singularity. Further from the element, however, the integrand is 
better approximated by a polynomial and the weights come closer to those of 
the standard rule. 
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Figure 2: Deviation of quadrature weights from corresponding Gauss-Legendre 
rule: AT = 64, M = 16. 

3.1 Accuracy 

The test of accuracy is to examine how well the new quadrature rule evaluates 
a reference integral. Table 1 shows the error in computing: 



for three values of R over the range 7r/64 < 6 < 3l7r/64. The case of 9 — tt/2 
was ignored because for odd n, I = which would not allow for a meaningful 
estimate of relative error which is defined: 



where K is the value of / estimated by numerical quadrature. 

Tables 1 and 2 show the error ei, incurred using the method of this paper, 
compared to £2 , the error using standard Gauss-Legendre quadrature, using low 
and high order rules. The calculation is carried out at three values of R to 
examine the change in error with distance from the clement and for various 
values of n. Table 1 shows the error for n = 0, . . . , 3, the important range for 
boundary element calculations, and a rule with A'' = 16 and M = 4. Prom the 
first two columns of errors, it is clear that the modified rule is far superior to 






(24) 
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Table 1: Root-mean-square error in reference integral computed with modified 
(ei) and standard (£2) quadratures at three values of R, N = 16, M = 4. 



R 




2-' 






2° 








2' 






n 


ei 




£2 


ei 




£2 




ei 




£2 







1.6 X 10" 


-12 


5.9 X 10" 


3.6 X 10" 


11 


1.4 X 10" 


-2 


2.1 X 10- 


-w 


4.5 X 10- 


-16 


1 


2.8 X 10" 


■13 


3.3 X 10° 


1.3 X 10" 


■10 


1.3 X 10- 


-2 


1.3 X 10- 


-16 


7.3 X 10- 


■16 


2 


6.0 X 10- 


-14 


1.9 X 10° 


1.0 X 10" 


-10 


1.3 X 10 


-2 


2.6 X 10- 


-16 


1.5 X 10- 


15 


3 


1.9 X 10" 


-13 


1.0 X 10° 


9.9 X 10- 


-11 


1.2 X 10" 


-2 


6.3 X 10- 


-16 


3.0 X 10- 


■15 



Table 2: Root-mean-square error in reference integral computed with modified 
(ei) and standard {€2} quadratures at three values of R, N = 64, M = 16. 



R 




2- 


■ 1 






2° 








2' 






n 


ei 




£2 




£1 




£2 




£1 




£2 







1.9 X 10- 


■10 


1.1 X 10" 


-01 


6.6 X 10- 


-15 


5.6 X 10- 


-12 


1.9 X 10" 


12 


3.7 X 10- 


■16 


3 


3.4 X 10- 


■10 


8.4 X 10" 


-03 


5.4 X 10- 


-15 


4.2 X 10- 


-12 


6.5 X 10- 


-13 


3.0 X 10- 


15 


6 


4.0 X 10- 


■10 


3.8 X 10" 


-03 


7.1 X 10- 


-15 


2.6 X 10- 


-12 


2.3 X 10- 


-12 


2.5 X 10- 


■14 


9 


4.0 X 10- 


■10 


8.0 X 10" 


-04 


8.5 X 10- 


-15 


9.9 X 10- 


-13 


1.0 X 10" 


12 


2.1 X 10- 


■13 


12 


4.0 X 10- 


■10 


1.4 X 10" 


-04 


8.4 X 10- 


-15 


6.6 X 10- 


-13 


2.8 X 10- 


-12 


1.7 X 10- 


12 


15 


3.8 X 10- 


■10 


2.2 X 10" 


-05 


8.3 X 10- 


-15 


2.3 X 10- 


-12 


5.2 X 10- 


-12 


1.3 X 10- 


■11 



the standard quadrature: its error is about twelve orders of magnitude lower 
than that of the Gaussian quadrature. Similarly, for i? = 1, the mean error is 
much smaller, being no worse than about 10"^°, rather than 10"^. Once the 
field point is far from the element, however, at i? = 2, both rules have similar 
accuracy. 

Table 2 shows similar results for a rule with = 64 and M = 16, compared 
to data for a standard Gauss-Legendre rule with A'^ = 64. Integrals of order up 
to n = 15 have been computed and, as before, at small R, the error behaviour of 
the new rule is orders of magnitude better than that of the standard technique. 
At larger R, however, the advantage is not so clear cut: at i? = 1, the error in 
the standard rule is around 10"^^, still larger than that from the method of this 
paper, but probably acceptable in many applications. When R = 2, the Gauss- 
Legendre rule gives results comparable to those of the new technique, although 
it performs better on low order polynomials. As might be expected, when the 
distance from the element is large enough, a high order Gauss-Legendre rule 
can capture enough of the behaviour of the integrand to accurately compute 
the integral. 

To examine the error behaviour in more detail, figures 3-5 show the relative 
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Present method R = 2 ^ 
Gauss-Lcgcndrc R = 2^^ 
Present method R = 2^ 
Gauss-Legendre R = 2^ 




7r/8 



7r/4 



37r/8 
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Figure 3: Relative error in logarithmically singular integral, A'' = 16, M = 4. 



error in computing: 



with error e defined as: 



where Jq is one of lQ°^\x,y), lQ'^'{x,y), lQ''{x,y) and Jq is the corresponding 
estimate using the quadrature rule. In each case, N = 1% and M = 4 and a 
sixteen point Gaussian quadrature was also applied for comparison. In each 
case, for large distances from the element, R = 2, both quadrature techniques 
are accurate, with errors of the order of machine precision. When R = 2~^, 
however, the error incurr(xi using Gaussian quadrature is unacccptably large, 
while the modified rule gives very accurate answers, again of the order of ma- 
chine precision. As might be expected, the error from the Gaussian quadrature 
is smaller as 9 7r/2, due to the greater distance from the clement, but it is 
never better than about 10~^, insufficient accuracy for most applications. 
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Present method R = 2 ^ 
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Figure 4: Relative error in near-singular integral, = 16, M = 4. 
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Figure 5: Relative error in near-hypersingular integral, A'' = 16, M = 4. 
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Figure 6: Relative error in near-hypersingular integral, A'' = 64, M = 16. 



Finally, figure 6 illustrates an interesting point about the error behaviour of 
the quadrature rule as the number of quadrature points is increased. It shows 
the error in the near-hypersingular integral /q {x,y) evaluated using a rule 
with TV = 64 and M = 16, with the error from a 64-point Gaussian quadrature 
shown for comparison. The first point is that the Gaussian quadrature is able 
to cope with the singularity for R = when 6 ^ n/8, because it can integrate 
polynomials of high enough order to be able to correctly handle the expansion 
of the integrand in Legendre polynomials. Also, as in the previous cases, it can 
correctly evaluate the integral for R = 2. 

The modified rule, however, has slightly worse error behaviour than for N = 
4, with the maximum error being higher for i? = and the error at J? = 2 
being greater across the full range of 0. The error is still small, being less than 
lO^^'', but the reason for the increase is imclear. It appears to be due to an 
ambiguity in expressing the integrand in terms of Legendre polynomials: the 
near-hypersingular part of the integrand f{t)/[{x — t)'^ + y^] is a ratio of two 
polynomials which can be written as a sum of a proper elementary function and 
a polynomial. This polynomial term is then represented twice in the quadrature 
rule, being handled by the unweighted Legendre polynomials, rrii, 1 < i < M 
in equation (14) and by the weighted polynomials. An interesting question for 
future developments of the technique will be how best to choose the elementary 
functions for the quadrature rule to give optimal accuracy for a given N. 
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4 Conclusions 



A method of deriving quadrature rules for the evaluation of the 'near-singular' 
integrals which arise in the boundary clement mcrthod has been derived. The 
performance of the technique has been assessed by evaluation of reference inte- 
grals and it has been found that it outperforms standard Gaussian quadrature 
rules with the same number of nodes for field points close to the clement. The 
error in the integral increases slightly with the number of points in the rule, a 
point which is to be investigated in future work. 

A Integrals of weighted polynomials 

To evaluate the required integrals of Legendre polynomials using the procedure 
of §2.2, we require a method of evaluating the integrals: 

I^^\x,y) =f_^^-^dt, (29a) 

4'^ i^, y) =1-1 [(.-t)Cy^]V. dt, (29b) 

7('°s)(a;,2/) = f',t^log[{x - + y^]^^ dt, (29c) 

Use of standard formulae [9, 2.171,2.263,2.728.1] yields: 

lL'H^,y) = '-^^+2xI^2i{^,y)-RH^^l,ix,y), (30a) 

/(i)(.x,y) = + ^:rli'\ix,y) ^^R^ ll^\{x,y), (30b) 

^aog)(^,y) = '"^^^+^17°^^- - ^J':\x,y) + ^I^:\x,y), (30c) 

where the recursions are seeded with: 

ll^\x,y) =i(tan-ii^-tan-il±3^), jf ^ (:r, y) = log |i + x/f (x, y), 

ll^\x,y) =log|^^, /W(x,2/)=i? + x/«(x,2/), 

and 

R = (a;2+y2-)i/2^ i?± = [(a;Tl)^ + y^]^/^ ^ = tan-iy/a;. 
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